Universality Class of Criticality in the Restricted Primitive Model Electrolyte 
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The 1:1 equisized hard-sphere electrolyte or restricted primitive model has been simulated via 
grand-canonical fine-discretization Monte Carlo. Newly devised unbiased finite-size extrapolation 
methods using temperature-density, (T,p), loci of inflections, Q = (m 2 ) 2 /(m 4 ) maxima, canonical 
and Cv criticality, yield estimates of (T c , p c ) to ±(0.04, 3)%. Extrapolated exponents and Q-ratio 
are (7, v, Q c ) = [1.24(3), 0.63(3); 0.624(2)] which support Ising (n = 1) behavior with (1.23 9 , 
O.63O3; 0.6236), but exclude classical, XY (n = 2), SAW (n = 0), and n=l criticality with potentials 
<fi(r) > $/r 4 ' 9 when r — > 00. 

PACS numbers: 02.70.Rr, 05.70.Jk, 64.60.Fr, 64.70.Fx 



Since the experiments of Singh and Pitzer in 1988 
||, an outstanding experimental and theoretical ques- 
tion has been: What is the universality class of Coulom- 
bic criticality? Early experimental data for electrolytes 
exhibiting phase separation driven by long-range ionic 
forces suggested classical or van der Waals (vdW) crit- 
ical behavior, with exponents /? = i, 7 = 1, 1/ = i, 
etc. jj], |^]: But the general theoretical consensus has been 
that asymptotic Ising-type criticality, with (3 ~ 0.326, 
7 ~ 1.239, v ~ 0.630, etc., should be expected @, |, §. 
Naively, one may argue that the exponential Debye 
screening of the direct ionic forces results in effective 
short-range attractions that can cause separation into 
two neutral phases: ion-rich and ion-poor |jl], ||, ||, ^, |j; 
the order parameter, namely, the ion density or concen- 
tration difference, is a scalar; so Ising-type behavior is 
indicated. Field-theoretic approaches support this pic- 
ture ||. 

However, the theoretical arguments are by no means 
rigorous and have not, so far, been tested by precise cal- 
culations for appropriate models. To do that is the aim 
of the researches reported here. We have studied a finely- 
discretized version |?]] of the simplest continuum model 
(considered by Debye and Huckel in 1923 |jl|, [2), three 
years before Ising's work), namely, the restricted prim- 
itive model (RPM), consisting of N = N + + A_ equi- 
sized hard spheres of diameter a, precisely half carrying 
a charge +qo and half — qo, in a medium (representing a 
solvent) of dielectric constant D. At a separation r > a, 
like (unlike) ions interact through the potential ztq^/Dr; 
thus appropriate reduced density, p = N/V for volume V, 
and temperature variables are 



p = pa 



T* = k B TDa/ql, t={T-T c )/T c . (1) 



Except at low densities and high temperatures, when 
the inverse Debye length n D a = (iirp* /T*) 1 / 2 is small, 
the RPM is intractable analytically or via series ex- 
pansions 0, ||, [|]. However, it has been much studied 
by Monte Carlo (MC) simulations fj], §, 0, |l) which 
have recently approached the consensus T* ~ 0.049, 



p* = 0.060-0.085. However, these values have been de- 
rived by assuming Ising-type criticality: on that basis 
Bruce-Wilding extrapolation procedures have been em- 
ployed |, (which, even then, neglect potentially im- 
portant, asymmetric 'pressure-mixing' terms |l2[].) It 
must be stressed that implementing appropriate finite- 
size extrapolation methods constitutes the heart of the 
computational task since a grand-canonical (GC) system 
confined in a simulation 'box' of dimensions L x L x L 
(with, say, periodic boundary conditions ]l3| ) cannot ex- 
hibit a sharp critical point; a finite canonical system may 
become critical but can display only classical or vdW be- 
havior Q. 

Thus, while previous RPM simulations ^ |l(| demon- 
strate consistency with Ising (or n = 1) behavior, no 
other universality classes are ruled out: see also 0, [lj, 
|l5| . Putative 'nearby' candidates are XY or n = 2 sys- 
tems (with 7 ~ 1.316, v ~ 0.670), self-avoiding walks 
(SAWs, n = 0: with 7 ~ 1.159, v ~ 0.588) @, || 
and long-range, \l r d +' T scalar systems (with d = 3, 
a < 2 — 77) 1 15, 17\. On the other hand, in a prepara- 
tory GCMC study II @(b)] of the hard -core square- 
well (HCSW) fluid — for which Ising criticality has long 
been anticipated — new, unbiased, finite-size extrapola- 
tion techniques enabled the n = 2 and classes to be 
convincingly excluded. 

Present approach. — We have now applied the meth- 
ods of II to the RPM; however, the extreme asymme- 
try of the critical region in the model (see Fig. [T]) has 
demanded further developments. By extending finite- 
size scaling theory |l8| and previous applications of the 
Binder parameter or fourth-moment ratio Jl7| , p"8[ [T9| 

Q L {T- p) = (m 2 ) 2 /(m 4 ) with m = p - (p) , (2) 

p0| to systems lacking symmetry, we have assembled ev- 
idence, outlined below, that excludes not only classical 
criticality in the RPM but also the XY and SAW univer- 
sality classes and (d = 3) long-range Ising criticality with 
a < 1.9. 

Our work employs multihistogram reweighting |plf and 
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FIG. 1: Approximate coexistence curve of the RPM in the 
(T, p) plane: open circles and fitted line. The estimated crit- 
ical point is shown as an uncertainty bar. The dashed curves 
are loci of CV(T) maxima at fixed p for L* = L/a = 8, 10, 
and 12. The loci labeled k = 1, and Q are explained in 
the text. The inset shows the canonical critical points T®{L"), 
p° c (L) (squares), and corresponding GC mean densities p c (L) 
(crosses) for L* = 9-12, the Cv{L) extrema T~(L),p~(L), 
for L* = 7-10 and 12 (solid circles), and the ^fp diameter, 
Py 2 (T), defined in the text (open squares). 



a (£ = 5)-level fine-discretization formulation (with a fine- 
lattice spacing a/C @)- Since £ < oo, nonuniversal pa- 
rameters, such as T*, will deviate slightly from their 
continuum limit (C — > oo) |t], ^2|; but, at this level, 
there are no serious grounds for contemplating changes 
in universality class. For the critical parameters we find 
T* = 0.05069(2) and p* = 0.0790(25): the confidence 
limits in parentheses refer, here and below, to the last 
decimal place quoted. The inset in Fig. [I] shows how these 
values are approached (i) by the canonical values T^(L), 
p°(L) and p\{L) (= (p)t°(l), m °(£) ©) derived from the 
isothermal density histograms [see II(2.18)-(2.23), Figs. 
1, 3], (ii) by T~(L) and p~{L), from the isochoric max- 
ima of CV(T; p; L) [see Fig. [j] and II Sec. Ill, Fig. 7], and 
(iii) by the yfp diameter, p 1 / 2 (T), denned below. 

Exponents 7 and v. — Before justifying the precision 
of our (T c , p c ) estimates, we consider their implications. 
The solid curves in Fig. || portray the effective susceptibil- 
ity exponent 7 ( J f (T 1 ; L) on the critical isochore above T c , 
as derived from xnn = V(m 2 ) = ksTp 2 Kx'- see 11(3.7). 
Within statistical precision the data are independent of 
the (T c , p c ) uncertainties. 

Also presented in Fig. |2| are the modified estimators 
7^jj(T) [defined as in 11(3.7) but with t replacing t'] eval- 
uated on the 'theta locus,' p#{T) = p c [tf+(l-tf)(T c /T)]. 
This relation approximates an effective symmetry lo- 
cus (II) above T c , derived from the behavior of the 
isothermal inflection loci pk(T;L), on which x^ = 
X NN(T,p;L)/p k is maximal [see II(2.26)-(2.32)]. The 
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FIG. 2: Effective susceptibility exponent 7^(T) for p — p c 
(solid curves) and %g(T) on the theta locus (dashed; see text), 
for sizes L* = 7-12 and 15. Values for vdW and for n = 0, 1, 
and 2 are marked on the 7 axis. 

k = 1 loci are shown in Fig. 1 for L* = L/a = 6, 8, 10, 12; 
the selected value $ — 0.20 corresponds roughly to k ~ 
0.60 (which may be identified with an optimal value: see 
II and [jl8|). However, the variation of the k loci when 
L increases is significantly more complicated in the RPM 
than in the HCSW fluid [11(c), 23]. 

Extrapolation of the effective susceptibility exponents 
in Fig. U and those on the k — locus, etc. [|TT|(c)1, to t = 
indicates 7 = 1.24(3), upholding Ising-type behavior 
while both XY and SAW values are implausible. 

To determine the exponent v we have examined the 
peak positions, Tj(L), of various properties, Yj (T: L), on 
the critical isochore. Finite-size scaling theory |18[ yields 
ATj (L) = T 3 (L) - T c ~ L~ x l v : Figure | demonstrates 
the estimation of \ jv (unbiased except for the imposed T c 
estimate) from the ratios ATj(Li)/ ATj(L,2) for various 
j (see |ll|(c)1), using an established approach [see 1(7)- 
(13), Fig. 1; 11(3.1)]. The data indicate v = 0.63(3), 
excluding classical but supportive of fsing (n = 1) criti- 
cality, while n = 2 and seem less probable. 

Estimation of T*. — Consider, now, Ql(T;p) in (§), 
when L — > 00. In any single-phase region of the (T, p) 
plane Ql — ► ^ , indicative of Gaussian fluctuations about 
(p); conversely, within a two-phase region, p-(T) < p < 
p_l_(T), one finds Ql — > 1 on the diameter, p(T) = ^(p_ + 
p+) for T < T C1 while, more generally, 

1 > Qoo(T; p) = 1 - 4y 2 /(l + 6y 2 + y 4 ) > \ , (3) 

where \y\ = 2|p - p(T)|/(p_ + p+) < 1. Finally, at 
criticality, Ql(T c ',Pc) approaches a universal value Q c 
which, for cubic boxes with periodic boundary condi- 
tions, is Q c = 0.4569 • • • for classical (vdW) @b)] or 
oo-range systems Qc)] but Q c (n=l) = 0.6236(2)^ for 
Ising fjj(d),(e)] and Q c {n = 2) = 0.8045(1) for XY ggf)] 
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FIG. 3: Estimation of the correlation exponent v from the 
deviations Tj (L) — T c for various properties Yj (T) on the crit- 
ical isochore: see text and (ll](c)]. Values for n — 0,1,2, 
i.e., SAW, Ising, and XY, and classical (vdW) criticality are 
indicated. 
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FIG. 4: Plots of Q L (T;p) on the Q-loci, pq(T;L), providing 
estimates for T c and Q c . Classical, XY and Ising values of Q 
are shown. 



systems, while Q c (n = 0) = |l9|(b)]. For long-range, 
l/r 3+<T systems, Q c (<t) and also 7(c), increase almost 
linearly from vdW to Ising values in the interval | < 
(7 < (j/u)„=i ~ 1.96 6 with Q c (cr = 1.9) ~ 0.600 and 
7 (a=1.9)~1.20 5 @b)]. 

The result (||) leads us to propose Q-loci, Pq{T; L), on 
which Ql(T; p) is maximal at fixed T. For T < T c these 
loci are observed to approach the diameter p(T) when L 
increases. (For T < T c , but not above T c , the Q-loci also 
follow the k = loci quite closely.) 

Figure ^ displays Ql(T; p) on the Q-loci Pq(T; L), for 
L* = 7-12. As often seen in plots for symmetric sys- 
tems Jl9| , inflection points and successive intersections, 
Tq(L), almost coincide! Scaling yields Ql{T c ;p c ) ~ 
L- e l v and |T c -T Q (i)| - with ip = [1 + 9)/ v, where 
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FIG. 5: Estimation of p* c from plots of (k — 0) and Q locus 
values at T c (open and solid squares) vs A/(L* + £o)^ for 
various values of ip and optimal shifts lo- The scale param- 
eter A has been invoked merely for graphical clarity. Note 
ip < 1.6 requires smaller shifts tending to exclude vdW criti- 
cality (ip = 2). 



(= aw) is the leading correction-to-scaling exponent; 
for classical and Ising criticality one has (6/l>, (p) = (1,3), 
~ (0.82, 2.41) (ie|. With this guidance, the large-scale in- 
set in Fig. | leads to our estimate T* ~ 0.05069(2) but 
also yields Q c ~ 0.624(2): this is surprisingly close to 
the Ising value Q and far from the vdW, XY, and SAW 
values — an unexpected bonus! Likewise, l/r 3+tT effective 
potentials with er < 1.9 are excluded. 

Estimation of p c . — Finally, we examine p%(L) and 
Pq(L), i.e., the (k = 0) and Q loci intersections with 
the estimated critical isotherm, T = T c . According to 
scaling, the deviations, Apg and Apg, decay as with 
ip = (X — a)/v Q, so we may suppose 1.2 < ip < 2 [jl6|. 
Figure^ displays the deviations vs for ip = 1.2, 1.4, 
1.7 and 2 with '4 shifts' [1(19), Fig. 2; 11(3.1)] chosen to 
provide linear plots. From these and further plots |ll](c)] 
we conclude p* = 0.0790(25). 

In further support of our p c estimate, we mention first 
that when the coexistence curve, p±(T), is plotted vs 
y/p* — as is reasonable since all powers p^ 2 for integral j 
appear in virial expansions for the RPM || — it becomes 
markedly more symmetrical [resembling (p, T) plots for 
the HCSW and other simple fluids]. Then, the corre- 
sponding diameter, \/[/o*^ 2 (T)] = \{\f~P*i- + \/p+)i is 
only mildly curved and naive extrapolations to T c yield 
p* = 0.078(4). 

In conclusion. — By implementing recently tested |14j 
and newly devised extrapolation techniques for nonsym- 
metric critical systems, our extensive grand-canonical 
Monte Carlo simulations for the RPM have provided, 
in toto, convincing evidence to exclude classical, XY 
(n = 2), or SAW (n = 0) critical behavior as well as 
long-range (effective) Ising interactions decaying more 
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slowly than l/r 4 ' 90 . Rather, the estimates for the expo- 
nents v and 7, and for the critical fourth- moment ratio, 
Q c , point to standard, short-range Ising-type criticality. 
Studies underway JTT](c)] should provide further confir- 
mation and additional quantitative results, such as the 
scale, Ro, of the equivalent single-component short-range 
attractions generated by the RPM near criticality. 
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of Energy, Office of Basic Energy Sciences (DE-FG02- 
01ER15121, AZP) is gratefully acknowledged. 
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